rm(list = ls())
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center',
                      fig.width = 8,
                      fig.height = 6,
                      dev = 'jpeg',
                      dpi = 300)
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo") 
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapido
library(here)
library(kableExtra)
library(ggthemes)
library(ggrepel)
library(ggridges)
#analisis
library(ggpubr)
library(easystats) # multiples unciones analiticas
library(sf)
library(tidyverse, quietly = TRUE)
library(modelsummary)
library(terra) # replace raster
library(TropFishR)
library(mixR)
library(CCAMLRGIS)

Background

Methodology

A important piece of information for a stock evaluation refers to the biological components such as average sizes and weights across areas and years. To do this, we will explore the biological data and prepare the output to add it into stock assessment integrate model (Methot & Wetzel, 2013).

The object ohbio2 come from data exploration analysis in data request CCAMLR data. This objetc have bio information from krill.

#cargo objeto
meta <- get("METADATA")
c1 <- get("C1")
ohbio <- get("OBS_HAUL_BIOLOGY")

Join data set with master as c1 set. This join is trought obs_haul_id variable to get geoposition variables

ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
names(ohbio2)
##  [1] "c1_id.x"                   "obs_haul_id"              
##  [3] "obs_logbook_id.x"          "obs_haul_number"          
##  [5] "haul_number.x"             "vessel_name"              
##  [7] "vessel_nationality_code"   "fishing_purpose_code"     
##  [9] "season_ccamlr"             "target_species"           
## [11] "asd_code"                  "trawl_technique"          
## [13] "catchperiod_code"          "date_catchperiod_start"   
## [15] "datetime_set_start"        "datetime_set_end"         
## [17] "datetime_haul_start"       "datetime_haul_end"        
## [19] "datetime_timezone"         "depth_gear_set_end_m"     
## [21] "depth_gear_haul_start_m"   "depth_bottom_set_end_m"   
## [23] "depth_bottom_haul_start_m" "latitude_set_end"         
## [25] "longitude_set_end"         "latitude_haul_start"      
## [27] "longitude_haul_start"      "gear_type_code"           
## [29] "gear_type"                 "mesh_code"                
## [31] "trawl_net_number"          "notes"                    
## [33] "trawl_duration_depth_h"    "trawl_duration_total_h"   
## [35] "krill_greenweight_kg"      "c1_id.y"                  
## [37] "obs_logbook_id.y"          "haul_number.y"            
## [39] "taxon_code"                "taxon_scientific_name"    
## [41] "taxon_family"              "maturity_stage"           
## [43] "sex_code"                  "length_total_cm"          
## [45] "greenweight_kg"

Firsts glance. Test how many register have by year. In this case, length_total_cm by season ccamlr. Same exercise in date period date_catchperiod_start to separate dates.

ohbio3 <- ohbio2 %>%
  mutate(Year = year(date_catchperiod_start),
         Month = month(date_catchperiod_start),
         Day = day(date_catchperiod_start)) %>% 
  #toupper() para convertir los valores a mayúsculas
  mutate(sex_code = toupper(sex_code))

Save data further analysis

length481 <-ohbio3 %>% 
  dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47, 43, 45) |>   filter(asd_code=="481")
#save(length481, file = "length481.RData")

First thing is get different rater layer to join krill data length according different porpoises.

# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast) 
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")
# Uso las agrupaciones de Strata
strata <- st_read("~/DOCAS/Mapas/Antarctic_SHPfiles/Strata.shp",
                quiet=T)
strata=st_transform(strata, "+proj=latlong +ellps=WGS84")

Test Strata SSMU, just to know another way to join data, but this kind of spatial structuration is deprecated to mamagemente use (Figure@ref(fig:maptest2).

Grouping bio data into stratas

names(length481)
##  [1] "vessel_nationality_code" "season_ccamlr"          
##  [3] "asd_code"                "trawl_technique"        
##  [5] "date_catchperiod_start"  "latitude_set_end"       
##  [7] "longitude_set_end"       "gear_type"              
##  [9] "maturity_stage"          "length_total_cm"        
## [11] "Year"                    "Month"                  
## [13] "sex_code"                "greenweight_kg"
ohbio6 <- st_as_sf(length481 %>% 
                     drop_na(latitude_set_end), 
                   coords = c("longitude_set_end", "latitude_set_end"),  
                  crs = "+proj=latlong +ellps=WGS84")

Study Area

# y testeo el mapa
ssmap <- ggplot()+
  geom_sf(data = strata |> 
            filter(ID != "Outer") |>  
           mutate(ID = str_replace_all(ID, "Extra", "GS")),color="red")+
  geom_sf(data = ohbio6 |>
            drop_na(length_total_cm) |> 
            filter(Year>2017), 
          aes(color = length_total_cm)) +
  geom_text_repel(data = strata |>  
            filter(ID != "Outer") |>  
           mutate(ID = str_replace_all(ID, "Extra", "GS")), 
            aes(x = Labx, y = Laby, 
                label = ID), 
            min.segment.length = 0,
                             box.padding = 2,
                             max.overlaps = 10)+
  geom_sf(data = coast2, colour="black", fill=NA)+
  scale_color_viridis_c(option = "G",
                        name="Length (cm)")+
  ylim(230000, 2220000)+
  xlim(-3095349 , -1858911)+
  # coord_sf(crs = 32610)+ #sistema de prpyecccion para campos completos
  coord_sf(crs = 6932)+
  theme_bw()
ssmap
Strata Maps in 48.1

Strata Maps in 48.1

Length composition by Strata CCAMLR to visualization first. First step is group data into to poligons strata.

strata <- st_make_valid(strata)
sf4 <- st_join(strata, ohbio6)
# Save an object to a file
saveRDS(sf4, file = "sf4.rds")

Load RData

# Restore the object
sf4 <- readRDS("sf4.rds")

Data

histogram length data to viz in anthor way.

jzstrata <- ggplot(sf4 %>% 
                     mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) %>% 
               filter(Year>2010,
                      ID !="Outer"),
             aes(x=length_total_cm, 
                 y = as.factor(Month), 
                 fill=ID))+
  geom_density_ridges(stat = "binline", bins = 30, 
                      scale = 1.9, 
                      draw_baseline = FALSE,
                      alpha=0.9)+
  facet_grid(Year~ID) +   
  geom_vline(xintercept = 3.6, color = "red")+
  scale_x_continuous(breaks = seq(from = 1, to = 10, 
                                  by = 2))+
  scale_y_discrete(breaks = seq(from = 1, 
                                to = 12, by = 4))+
  scale_fill_viridis_d(name="Strata",
                       option="F")+
  theme_few()+
  xlab("Length (cm.)")+
  ylab("")
jzstrata

Now, must filter the DF

sf5 <- sf4 |> 
  dplyr::select(c("ID",
        "date_catchperiod_start",
        "length_total_cm",
        "Year",
        "Month",
        "sex_code",
         "greenweight_kg")) |> 
  mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) |>  
  filter(ID !="Outer") |> 
  data.frame()

Calculating the proportion of records per year and per stratum in a table. This table allows us to establish criteria for defining the strata on which we will analyze krill parameter. We leave out just JOIN.

# Calcular la proporción
proporcion <- (round(table(sf5$ID) / sum(table(sf5$ID))*100,2))

# Crear una tabla con las proporciones
tabla_proporcion <- as.data.frame(proporcion)
tabla_proporcion$Categoria <- rownames(tabla_proporcion)
rownames(tabla_proporcion) <- NULL

# Renombrar las columnas
colnames(tabla_proporcion) <- c("Strata", "%")

# Mostrar la tabla utilizando kbl()
kbl(tabla_proporcion, 
    caption = "Proportion length register by Strata")  |> 
  kable_classic(full_width = F, 
                html_font = "Cambria") |> 
  kable_styling(bootstrap_options = "striped", 
                latex_options = "striped")
Proportion length register by Strata
Strata % NA
BS 62.71 1
EI 6.58 2
GERLASHE 12.13 3
JOIN 0.14 4
SSIW 18.43 5

Parameters estimation k and L_inf_

Method ELEFAN

sf5fil <- sf5 |>
  filter(ID %in% c("BS", "EI", "GERLASHE", "SSIW")) |>  
  mutate(date_catchperiod_start = as.Date(date_catchperiod_start)) |>  
  mutate(yearly_Group = floor_date(date_catchperiod_start, "year")) |>   # New column
  drop_na(length_total_cm) 

# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")
# Crear una lista para almacenar los resultados de cada ID
lfq_results <- list()
# Iterar sobre cada ID
for (id in unique(sf5fil$ID)) {
  # Obtener el nombre correspondiente
  nombre <- nombres[which(unique(sf5fil$ID) == id)]
  # Filtrar los datos para el ID actual
  df_id <- sf5fil %>% filter(ID == id)
  # Crear un objeto lfq para el ID actual
  lfq <- lfqCreate(data = df_id,
                   Lname = "length_total_cm",
                   Dname = "yearly_Group",
                   bin_size = 0.1)
  # Agregar el resultado a la lista
  lfq_results[[id]] <- lfq
  # Graficar el objeto lfq
  plot(lfq, Fname = "catch",
       main = nombre)
}

Identified distribution

ei <- Bhattacharya(lfq_results$EI)
## Interactive session needed for Bhattacharya.
ss <- Bhattacharya(lfq_results$SSIW)
## Interactive session needed for Bhattacharya.
bs <- Bhattacharya(lfq_results$BS)
## Interactive session needed for Bhattacharya.
gs <- Bhattacharya(lfq_results$GERLASHE)
## Interactive session needed for Bhattacharya.

ahora damos asignaciones a curvas para representación

# Crear una lista para almacenar los resultados de PW_results <- list()

# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")

# Iterar sobre cada objeto lfq almacenado en lfq_results
for (i in seq_along(lfq_results)) {
  lfq <- lfq_results[[i]]
  
  # Obtener el nombre correspondiente
  nombre <- nombres[i]
  
  # Restructurar el objeto lfq
  lfqbin <- lfqRestructure(lfq, MA = 3, addl.sqrt = TRUE)
  
  # Graficar el objeto lfq reestructurado
  plot(lfqbin, hist.col = c("white", "black"),
       image.col = c(rep(rgb(1,0.8,0.8),1000), "white", 
                     rep(rgb(0.8,0.8,1),1000)),
       ylim = c(0,max(lfqbin$midLengths+0.5)),
       main = nombre)
  
  # Ajustar curvas al objeto lfq reestructurado
  tmp <- lfqFitCurves(lfqbin, par = list(Linf=6.5, 
                                          K=0.45,
                                          t_anchor=0.5),
                      draw = TRUE, col=4, lty=2)
}

Proceed to fit Modal Progression analysis with ELEFAN

Parameters to SSIW

# Response surface analyss
res_RSAsswi <- ELEFAN(lfq = lfq_results$SSIW,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$SSIW$dates[3],
                  cross.midLength = lfq_results$SSIW$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAsswi <- ELEFAN_SA(lfq_results$SSIW, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.005021 maxTime = 30
## Emini is: -0.2844106025
## xmini are:
## 5.324971939 0.912516579 0.6134328698 0.09638525543 0.3336850045 
## Totally it used 30.005036 secs
## No. of function call is: 1649
# run ELEFAN with genetic algorithm
res_GAsswi <- ELEFAN_GA(lfq_results$SSIW, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

Table with different method

# show results
RSAsswi <- unlist(res_RSAsswi$par)
GAsswi <- unlist(res_GAsswi$par)
SAsswi <- unlist(res_SAsswi$par)
parsswi <- round(rbind(RSAsswi, GAsswi, SAsswi),3)

parsswi %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH to SSWI") %>%
    kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE)
Parametres LH to SSWI
Linf K t_anchor C ts phiL
RSAsswi 5.448 0.640 0.790 0.500 0.500 1.279
GAsswi 5.989 0.579 0.532 0.859 0.880 1.318
SAsswi 5.325 0.913 0.613 0.096 0.334 1.413

Parameters to EI

# Response surface analyss
res_RSAei <- ELEFAN(lfq = lfq_results$EI,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$EI$dates[3],
                  cross.midLength = lfq_results$EI$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAei <- ELEFAN_SA(lfq_results$EI, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.001442 maxTime = 30
## Emini is: -0.2111109108
## xmini are:
## 5.863710068 0.9813904597 0.873453021 0.9993372839 0.5819024704 
## Totally it used 30.001453 secs
## No. of function call is: 2048
# run ELEFAN with genetic algorithm
res_GAei <- ELEFAN_GA(lfq_results$EI, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAei <- unlist(res_RSAei$par)
GAei <- unlist(res_GAei$par)
SAei <- unlist(res_SAei$par)
parei <- round(rbind(RSAei, GAei, SAei),3)

Parameters to BS

# Response surface analyss
res_RSAbs <- ELEFAN(lfq = lfq_results$BS,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$BS$dates[3],
                  cross.midLength = lfq_results$BS$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAbs <- ELEFAN_SA(lfq_results$BS, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.00739 maxTime = 30
## Emini is: -0.3184663901
## xmini are:
## 5.310221989 0.8441457965 0.2786562045 0.3598986657 0.9663543715 
## Totally it used 30.007403 secs
## No. of function call is: 1829
# run ELEFAN with genetic algorithm
res_GAbs <- ELEFAN_GA(lfq_results$BS, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAbs <- unlist(res_RSAbs$par)
GAbs <- unlist(res_GAbs$par)
SAbs <- unlist(res_SAbs$par)
parbs <- round(rbind(RSAbs, GAbs, SAbs),3)

Parameters to GERLASHE

# Response surface analyss
res_RSAgs <- ELEFAN(lfq = lfq_results$GERLASHE,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$GERLASHE$dates[3],
                  cross.midLength = lfq_results$GERLASHE$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAgs <- ELEFAN_SA(lfq_results$GERLASHE, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.004488 maxTime = 30
## Emini is: -0.3892859855
## xmini are:
## 5.28623276 0.9116002417 0.3847553197 0.06077500898 0.3120612055 
## Totally it used 30.0045 secs
## No. of function call is: 2521
# run ELEFAN with genetic algorithm
res_GAgs <- ELEFAN_GA(lfq_results$GERLASHE, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAgs <- unlist(res_RSAgs$par)
GAgs <- unlist(res_GAgs$par)
SAgs <- unlist(res_SAgs$par)
pargs <- round(rbind(RSAgs, GAgs, SAgs),3)

total Table with L_inf_ and K

t_k_linf <- rbind(parsswi, parei, parbs, pargs)
t_k_linf %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH in different strata with 3 types of algoritms") %>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE)
Parametres LH in different strata with 3 types of algoritms
Linf K t_anchor C ts phiL
RSAsswi 5.448 0.640 0.790 0.500 0.500 1.279
GAsswi 5.989 0.579 0.532 0.859 0.880 1.318
SAsswi 5.325 0.913 0.613 0.096 0.334 1.413
RSAei 5.879 0.470 0.750 0.500 0.500 1.211
GAei 5.301 0.891 0.435 0.359 0.570 1.399
SAei 5.864 0.981 0.873 0.999 0.582 1.528
RSAbs 5.448 0.640 0.790 0.500 0.500 1.279
GAbs 6.032 0.673 0.420 0.431 0.811 1.389
SAbs 5.310 0.844 0.279 0.360 0.966 1.377
RSAgs 5.448 0.640 0.790 0.500 0.500 1.279
GAgs 5.518 0.593 0.611 0.732 0.454 1.256
SAgs 5.286 0.912 0.385 0.061 0.312 1.406

Method NO LINEAR MODEL

We use firts a descomposition methot with library mixR and then we use nlmer (Bates et al. (2015)) to calcularte parametrs

# Agregar la nueva columna con los intervalos de longitud total
sf6 <- sf5 |> 
  drop_na(length_total_cm) |> 
  filter(ID %in% "GERLASHE",
         Year==2018,
         Month==5)

sf5a <-  sf5 |> 
  drop_na(length_total_cm ) |> 
  filter(ID %in% "BS",
         Year==2018)
  
# compruebo el n de comp

s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')

s_normal = select(sf5a$length_total_cm, ncomp = 2:6)




# fit a Normal mixture model
mod1 = mixfit(sf5a$length_total_cm, ncomp = 4,
              pi = c(0.5, 0.5, 0.5, 0.5), 
              mu = c( 3, 4, 5, 6), 
              sd = c(0.2, 0.2, 0.2, 0.2))
mod2_weibull = mixfit(sf6$length_total_cm, family = 'weibull', ncomp = 3)

s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')

s_normal = select(sf5a$length_total_cm, ncomp = 2:6)
plot(s_weibull)
plot(s_normal)

# plot the fitted model# plot the fitted model# plot the fitted model
plot(mod1)
plot(mod2_weibull)

# fit a Normal mixture model (equal variance)
mod1_ev = mixfit(sf6$length_total_cm, ncomp = 2, ev = TRUE)

Natural Mortality in Krill

In this exercises different bioanalogic methods are tested with General Algoritm

MSSWI <- c(res_GAsswi$par, list(agemax = res_RSAsswi$agemax))
# use the function M_empirical to estimate natural mortality
Msswi <- M_empirical(Linf = MSSWI$Linf, K_l = MSSWI$K, 
                  tmax = MSSWI$agemax, temp = 5,
                  tm50 = 2,
                  method = c("Pauly_Linf", "Hoenig", 
                            "Hoenig",
                             "AlversonCarney",
                             "RikhterEfanov"
                             ))
# Bransfield
BS <- c(res_GAbs$par, list(agemax = res_RSAbs$agemax))

# use the function M_empirical to estimate natural mortality
Mbs <- M_empirical(Linf = BS$Linf, K_l = BS$K, 
                  tmax = BS$agemax, temp = 5,
                  tm50 = 2,
                  method = c("Pauly_Linf", "Hoenig", 
                              "Hoenig",
                             "AlversonCarney",
                             "RikhterEfanov"
                             ))
# GERLASHE
GS <- c(res_GAgs$par, list(agemax = res_RSAgs$agemax))

# use the function M_empirical to estimate natural mortality
Mgs <- M_empirical(Linf = GS$Linf, K_l = GS$K, 
                  tmax = GS$agemax, temp = 5,
                   tm50 = 2,
                  method = c("Pauly_Linf", "Hoenig", 
                             "Hoenig",
                             "AlversonCarney",
                             "RikhterEfanov"
                             ))
# Elephand Island
MEI <- c(res_GAei$par, list(agemax = res_RSAei$agemax))

# use the function M_empirical to estimate natural mortality
Mei <- M_empirical(Linf = MEI$Linf, K_l = MEI$K, 
                  tmax = MEI$agemax, temp = 5,
                  tm50 = 2,
                  method = c("Pauly_Linf", "Hoenig", 
                             "Hoenig",
                             "AlversonCarney",
                             "RikhterEfanov"
                             ))
# junto las bases

Total_M <- cbind(Mei, Mbs, Mgs, Msswi)
colnames(Total_M) <- c("EI", "BS" , "GS", "SSWI")
mean <- colMeans(Total_M, na.rm = TRUE)
# Agregar los promedios como una nueva fila al final del array
Total_M_Mean <- rbind(Total_M, Mean = mean)
Total_M_Mean  %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Estimated M by Strata") %>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE)
Estimated M by Strata
EI BS GS SSWI
Alverson and Carney (1975) 0.6030 0.7790 0.8530 0.8660
Hoenig (1983) - Joint Equation 0.8690 0.8690 0.8690 0.8690
Hoenig (1983) - Fish Equation 0.8470 0.8470 0.8470 0.8470
Pauly (1980) - Length Equation 1.2090 0.9710 0.9160 0.8820
Rikhter and Efanov (1976) 0.7680 0.7680 0.7680 0.7680
NA NA NA NA NA
Mean 0.8592 0.8468 0.8506 0.8464

References

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Methot, R. D., & Wetzel, C. R. (2013). Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research, 142, 86–99. https://doi.org/10.1016/j.fishres.2012.10.012